

# Replication files for JOP article:
# Beyond the glass ceiling, more `housework'? 
# Womens' work assignment, performance and influence in political institutions


# By Louisa Boulaziz 


# Description of script: This script includes all the material needed to reproduce
# the results presented in the manuscript (the article).

#----------------------------------------------------------------------
# PACKAGES

library(tidyverse)
library(gridExtra)
library(stargazer)
library(mclogit)
library(survival)
library(survminer)
library(coxed)
library(eha)
library(MASS)
library(lmtest)
library(multiwayvcov)
library(marginaleffects)

#----------------------------------------------------------------------
# REPLICATING RESULTS FROM TABLE 2 IN THE ARTICLE 

# Loading in the choice-set data

load("data-cs.rdata")

# Estimating model 1

mod_ms <- clogit(is_judge_rapporteur ~ 
                   # gender
                   is_female + 
                   # judge level at case variables 
                   years_in_court_when_case_decided+
                   +age_judge_when_case_decided+ decisions_per_year+
                   # background
                   former_gc_judge + former_ag + was_judge + was_academic +
                   was_civil_servant 
                 + was_lawyer + was_politician + member_state
                 # strata
                 + strata(iuropa_decision_id),
                 # data
                 data = data)

# Estimating model 2

mod_salience_ms <- mclogit(cbind(is_judge_rapporteur, set) ~ 
                             
                             # gender
                             is_female* salience+
                             
                             # judge level at case variables 
                             years_in_court_when_case_decided + age_judge_when_case_decided +
                             decisions_per_year+
                             # background
                             former_gc_judge + former_ag + was_judge + was_academic +
                             was_civil_servant + was_lawyer + was_politician + member_state,
                           
                           data = data)

# Producing results table for article 

stargazer(mod_ms, mod_salience_ms, type = "latex", omit= c("member_state"), single.row = T,
          dep.var.caption= "Likelihood of case assignment",
          dep.var.labels = "",
          omit.stat = c("rsq", "ll", "logrank", "wald", "max.rsq", "lr"),
          covariate.labels = c("Female judge", 
                               "Case importance", 
                               "Years in court", 
                               "Age", 
                               "Work load", 
                               "Former GC judge", 
                               "Former AG", 
                               "Former judge", 
                               "Former academic", 
                               "Former civil servant", 
                               "Former lawyer", 
                               "Former politician", 
                               "Female x Case importance"
          ),
          add.lines=list(c("Member state fixed effects","YES", "YES")), 
          header = FALSE, 
          font.size = "tiny", 
          title = "Effect of gender and case importance on case assignments in the CJEU")


#----------------------------------------------------------------------
# REPLICATING RESULTS FROM TABLE 1 IN THE ARTICLE 

# Loading in the data

load("time_data4.Rdata")

time_data2 <- ddd

# Model 1
duration00 <- glm.nb(duration_days ~
                       is_female+
                       # MS fixed effects 
                       member_state,
                     data = time_data2)


duration00$ses  <- coeftest(duration00,
                            vcov. = cluster.vcov(duration00,
                            cluster = time_data2$iuropa_decision_id))[,2]


duration00$ps  <- coeftest(duration00,
                           vcov. = cluster.vcov(duration00,
                           cluster = time_data2$decision_year))[,4]

# Model 2
duration0 <- glm.nb(duration_days ~
                      is_female * salience +
                      # MS fixed effects 
                      member_state,
                    data = time_data2)


duration0$ses  <- coeftest(duration0,
                           vcov. = cluster.vcov(duration0,
                           cluster = time_data2$iuropa_decision_id))[,2]


duration0$ps  <- coeftest(duration0,
                          vcov. = cluster.vcov(duration0,
                          cluster = time_data2$decision_year))[,4]

# Model 3
duration1 <- glm.nb(duration_days ~
                      is_female + salience + 
                      # Background 
                      former_gc_judge + former_ag + was_judge + was_academic + 
                      was_civil_servant + was_lawyer + was_politician + 
                      # Judge experience in court and age 
                      age_judge_when_case_decided + years_in_court+ decisions_per_year + 
                      # Case complexity 
                      count_policy_areas + count_cited_documents + hearing +
                      decision_type +  test_policy_area + 
                      # MS fixed effects 
                      member_state,
                    data = time_data2)

duration1$ses  <- coeftest(duration1,
                           vcov. = cluster.vcov(duration1,
                           cluster = time_data2$iuropa_decision_id))[,2]


duration1$ps  <- coeftest(duration1,
                          vcov. = cluster.vcov(duration1,
                          cluster = time_data2$decision_year))[,4]


# Model 4

duration2 <- glm.nb(duration_days ~
                      is_female * salience + 
                      # Background 
                      former_gc_judge + former_ag + was_judge + was_academic + 
                      was_civil_servant + was_lawyer + was_politician + 
                      # Judge experience in court and age 
                      age_judge_when_case_decided + years_in_court+ decisions_per_year + 
                      # Case complexity 
                      count_policy_areas + count_cited_documents + hearing +
                      # Type of case 
                      decision_type +  test_policy_area+ 
                      # MS fixed effects 
                      member_state,
                    data = time_data2)



duration2$ses  <- coeftest(duration2,
                           vcov. = cluster.vcov(duration2,
                           cluster = time_data2$iuropa_decision_id))[,2]


duration2$ps  <- coeftest(duration2,
                          vcov. = cluster.vcov(duration2,
                          cluster = time_data2$decision_year))[,4]

# Producing table 

stargazer(duration0, duration00, duration1, duration2, type = "latex", omit = c("member_state", "test_policy_area"), dep.var.caption= "Duration days",
          dep.var.labels = "",
          single.row =T, 
          omit.stat = c("rsq", "ll", "logrank", "wald", "max.rsq", "lr"),
          covariate.labels = c("Female judge", 
                               "Case importance", 
                               "Case with oral hearing",
                               "Former GC judge", 
                               "Former AG", 
                               "Former judge", 
                               "Former academic", 
                               "Former civil servant", 
                               "Former lawyer", 
                               "Former politician", 
                               "Age", 
                               "Years in court", 
                               "Workload",
                               "Count policy areas",
                               "Count cited documents", 
                               "Case with hearing",
                               "Procedure, Preliminary ruling",
                               "Female judge x Case importance"),
          add.lines=list(c("Member state fixed effects",
                           "Yes","Yes", "Yes", "Yes"), 
                         c("Policy areas dummies","No", "No", "Yes", "Yes")), 
          header = FALSE, 
          font.size = "tiny", 
          title = "Effect of gender on case duration, with and without member state fixed effects")




#----------------------------------------------------------------------
# PRODUCING FIGURE 3 FROM MODELS AS SHOWN IN THE ARTICLE 

# Model 3 

p <- avg_predictions(duration1, by = "is_female")

model_effects1 <- ggplot(p, aes(x = as.factor(is_female), y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_pointrange() + 
  coord_flip() + 
  scale_x_discrete(labels=c("1" = "Female judge", "0" = "Male judge"))+
  xlab("") + 
  ylab("Case duration in days") + 
  ggtitle("Gender difference in case duration for all cases") + 
  theme_classic()

model_effects1

## Model 4 

p1 <- avg_predictions(duration2, by = c("is_female", "salience"))


p1 <- p1 %>% 
  filter(salience == 1) 

model_effects2 <- ggplot(p1, aes(x = as.factor(is_female), y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_pointrange() + 
  coord_flip()+
  scale_x_discrete(labels=c("1" = "Female judge", "0" = "Male judge"))+
  xlab("") + 
  ylab("Case duration in days") + 
  ggtitle("Gender difference in case duration for important cases only") + 
  theme_classic()
model_effects2

#----------------------------------------------------------------------
# REPRODUCING FIGURE 1 AND 2 FROM THE ARTICLE 

# Female Chief justices on high courts


load("Percentage_of_women_chief_justices_on_high_courts.Rdata")

data1 <- data1 %>%
  filter(country == "European Union - 27 countries (from 2020)")


plot2 <- data1 %>% 
  filter(country == "European Union - 27 countries (from 2020)") %>%
  rename(Area = country) %>%
  ggplot(aes(as.numeric(year), percentage_women, group= Area)) + 
  geom_path(aes(linetype = Area))+
  ggtitle("Percentage of women chief justices in national high courts")+ 
  xlab("Year") + 
  ylab("Percentage") + 
  theme_classic()

plot2

# Women in national high courts

load("Percentage_of_women_on_high_courts.Rdata")


plot1 <- data1 %>% 
  filter(year == 2023 & country != "European Union - 28 countries (1993-2020)") %>%
  ggplot(aes(reorder(country, -percentage_women), percentage_women)) + 
  geom_col() + 
  geom_hline(yintercept = 50, linetype = "dashed", color = "red")+
  coord_flip() + 
  #ggtitle("Percentage of women judges in national 
  # high courts in 2023")+ 
  xlab("Country/Area") + 
  ylab("Percentage") + 
  theme_classic()

plot1

# Women in international courts

load("Percentage_of_men_and_women_on_ICs.Rdata")

complete$Time <- ifelse(complete$Time == "Court of Justice of the European Union", "x", complete$Time)

complete$Time <- ifelse(complete$Time == "General Court", "General Court of the European Union", complete$Time)

complete$Time <- ifelse(complete$Time == "European Court of Justice", "Court of Justice of the European Union", complete$Time)

complete <- complete %>% 
  rename(Court = Time)

plot3 <- complete %>% 
  filter(Court != "x") %>% 
  ggplot(aes(as.numeric(year), as.integer(percentage_women), group= Court, color = Court)) + 
  geom_path(aes(linetype = Court))+
  ggtitle("Percentage of women judges on international courts")+ 
  xlab("Year") + 
  ylab("Percentage") + 
  theme_classic() 

plot3

